{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regularized Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn import datasets\n",
    "\n",
    "boston = datasets.load_boston()\n",
    "X = boston['data']\n",
    "y = boston['target']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before building the `RegularizedRegression` class, let's define a few helper functions. The first function standardizes the data by removing the mean and dividing by the standard deviation. This is the equivalent of the `StandardScaler` from `scikit-learn`.\n",
    "\n",
    "The `sign` function simply returns the sign of each element in an array. This is useful for calculating the gradient in Lasso regression. The `first_element_zero` option makes the function return a 0 (rather than a -1 or 1) for the first element. As discussed in the {doc}`concept section </content/c2/s1/regularized>`, this prevents Lasso regression from penalizing the magnitude of the intercept. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def standard_scaler(X):\n",
    "    means = X.mean(0)\n",
    "    stds = X.std(0)\n",
    "    return (X - means)/stds\n",
    "\n",
    "def sign(x, first_element_zero = False):\n",
    "    signs = (-1)**(x < 0)\n",
    "    if first_element_zero:\n",
    "        signs[0] = 0\n",
    "    return signs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `RegularizedRegression` class below contains methods for fitting Ridge and Lasso regression. The first method, `record_info`, handles standardization, adds an intercept to the predictors, and records the necessary values. The second, `fit_ridge`, fits Ridge regression using\n",
    "\n",
    "$$\n",
    "\\hat{\\bbeta} = \\left( \\bX^\\top \\bX + \\lambda I'\\right) ^{-1} \\bX^\\top \\by.\n",
    "$$\n",
    "\n",
    "The third method, `fit_lasso`, estimates the regression parameters using gradient descent. The gradient is the derivative of the Lasso loss function:\n",
    "\n",
    "$$\n",
    "\\dadb{L(\\bbetahat)}{\\bbetahat} = - \\bX^\\top\\left( \\by - \\bX \\bbetahat \\right) + \\lambda I'\\text{ sign}(\\bbetahat).\n",
    "$$\n",
    "\n",
    "The gradient descent used here simply adjusts the parameters a fixed number of times (determined by `n_iters`). There many more efficient ways to implement gradient descent, though we use a simple implementation here to keep focus on Lasso regression. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "class RegularizedRegression:\n",
    "        \n",
    "    def _record_info(self, X, y, lam, intercept, standardize):\n",
    "        \n",
    "        # standardize \n",
    "        if standardize == True: \n",
    "            X = standard_scaler(X)\n",
    "        \n",
    "        # add intercept\n",
    "        if intercept == False: \n",
    "            ones = np.ones(len(X)).reshape(len(X), 1) # column of ones \n",
    "            X = np.concatenate((ones, X), axis = 1) # concatenate\n",
    "            \n",
    "        # record values\n",
    "        self.X = np.array(X)\n",
    "        self.y = np.array(y)\n",
    "        self.N, self.D = self.X.shape\n",
    "        self.lam = lam\n",
    "        \n",
    "    def fit_ridge(self, X, y, lam = 0, intercept = False, standardize = True):\n",
    "        \n",
    "        # record data and dimensions\n",
    "        self._record_info(X, y, lam, intercept, standardize)\n",
    "        \n",
    "        # estimate parameters\n",
    "        XtX = np.dot(self.X.T, self.X)\n",
    "        I_prime = np.eye(self.D)\n",
    "        I_prime[0,0] = 0 \n",
    "        XtX_plus_lam_inverse = np.linalg.inv(XtX + self.lam*I_prime)\n",
    "        Xty = np.dot(self.X.T, self.y)\n",
    "        self.beta_hats = np.dot(XtX_plus_lam_inverse, Xty)\n",
    "        \n",
    "        # get fitted values\n",
    "        self.y_hat = np.dot(self.X, self.beta_hats)\n",
    "        \n",
    "        \n",
    "    def fit_lasso(self, X, y, lam = 0, n_iters = 2000,\n",
    "                  lr = 0.0001, intercept = False, standardize = True):\n",
    "\n",
    "        # record data and dimensions\n",
    "        self._record_info(X, y, lam, intercept, standardize)\n",
    "        \n",
    "        # estimate parameters\n",
    "        beta_hats = np.random.randn(self.D)\n",
    "        for i in range(n_iters):\n",
    "            dL_dbeta = -self.X.T @ (self.y - (self.X @ beta_hats)) + self.lam*sign(beta_hats, True)\n",
    "            beta_hats -= lr*dL_dbeta \n",
    "        self.beta_hats = beta_hats\n",
    "        \n",
    "        # get fitted values\n",
    "        self.y_hat = np.dot(self.X, self.beta_hats)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell runs Ridge and Lasso regression for the {doc}`Boston housing</content/appendix/data>` dataset. For simplicity, we somewhat arbitrarily choose $\\lambda = 10$—in practice, this value should be chosen through cross validation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set lambda\n",
    "lam = 10\n",
    "\n",
    "# fit ridge\n",
    "ridge_model = RegularizedRegression()\n",
    "ridge_model.fit_ridge(X, y, lam)\n",
    "\n",
    "# fit lasso\n",
    "lasso_model = RegularizedRegression()\n",
    "lasso_model.fit_lasso(X, y, lam)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below graphic shows the coefficient estimates using Ridge and Lasso regression with a changing value of $\\lambda$. Note that $\\lambda = 0$ is identical to ordinary linear regression. As expected, the magnitude of the coefficient estimates decreases as $\\lambda$ increases. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDoAAAJ1CAYAAAA485A6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdebgjZZn38e/d3ewIzabiMjYoCgIziuIyg9gCihuDguPoINAuOOCrDiqiiAq4IzAKAw7CqC2Lu4CioKjQCuICIgKy7yogm81qg8Dz/vHUsdPpnCSVk5yq1Pl+rquu06k8qbpTqfzOue6uJVJKSJIkSZIkNcGsqguQJEmSJEkaFhsdkiRJkiSpMWx0SJIkSZKkxrDRIUmSJEmSGsNGhyRJkiRJagwbHZIkSZIkqTFsdEiSJEmSpMaw0SFJkiRJkhrDRkcfImJ+RKSIWNDn+OsjYtFoqxofEbFBRJwSEbcV23Fht+fKbu+2dQ382iYYYF9dUIyfP4R1nxsRD0TELyNi3lSXJ02FuT015vb0MbelzNyeGnN7+pjb42FGNjpads7W6d6IuCAi3hURc6qucVgiYtWI2Dsizo6IOyPibxHx54g4rfjSTcd7XQi8EDgY2BX4fJ/P1VpEPCMiDqx7yExjnYcBxwHPBfYZ8br6FhH7RcQ3I+La4rt+fZexs4oMuDwilkTEHyLisIhYbarjR7nsmcDcNreHwdxejrk9w7N1lMxtc3sYzO3l1DW365etKaUZNwHzgQR8BXgD+Qv/XuDiYv4xbeNnASsDs/tc/vXAohq8z6cAVxTv6UfFe3wj8L7icQI+PeIaVgIeAY7o97my23tYrx1gXQuKbTi/6s+62/vvVuew3wMwB7gH+EXV26KlpgTcUezzdwLXdxl7eDH+JGAP4L+BvwFnArOmMn6Uy54Jk7ltbg/pvZnbyy/P3J7B2Triz9HcNreH8d7M7eWXV8fcrl22NqaTOqALUkonTDyIiM8BlwNviYj9U0q3AaSUHgGWVFTjQCJiFeB7wIbAzimlk9qGHBwRWwJbjriUxwBB/kOlr+emsr3H8bMapqrff0rpoYi4BNgsIiIVyVexJ6eUrgUoalu906CI2BR4B3BSSmnnlvnXAUcAryP/sVZ6/CiXPQOZ2+Z2o1T9/s1ts3UamNvmdqNU/f7rltu1zdaquz9VTCztMO/T4blvFs89r8P4BW1jnwh8A7gLuBs4FXgyHTrMwDzg28W4u4DvABt0GluMXwn4APB78hdpcbH8Z/b5Ht9R1PypEttlXeAo4A/Ag8XPo4B1BqmPfJhc6jDN7/HcZNt7RWBf4ELg/mI7ng+8vY/Pqq/tydKO6zbkw8GuAR4ArgR2bxl34CT1LyyeX7kYc0VR62Ly/2Ac0sfn8KRiWQe2zT+jmL932/xfAZd2ev991NnX+y2xDwXwm2KZG1T9Xe9Q3yVM8j+DwMeKul/QNn9l4D7gtEHHj3LZM2XC3J7sNea2uW1um621nDC3J3uNuW1uNyq3qWm2zvQjOjp5cvGzU0f07yJiLvAzcvgeDVxKPvftLGCVtrHrAGeTO6pHA5cBLyjGdjqXfwXgB8A/A8cDRwJrkg8D+nlEbJ1SOr/H+3hN8fOYHuMm1rkmcC758LsvAhcAzwT2AraJiOeklO4pWd/nySH5GeBk8qFMFO+/23ObdKhvReCH5GA5AziBHKCbAzsVNUz23gbZnp8gf46fJwfRXsDCiLg6pfTzot71gbcWYy8rXndN8fMo4E3kc+g+A8wGNiIHXFcppRuKDui25OCceP//Qj70cFvgs8X8NYBnkferTnrV2e/77ddewBbFvzcHrus2OCJmAWuXWP6dKXfRR2FL8vb9devMlNKSiLiQ5f83psz4US5b5ra5nZnb5jZgto4Jc9vcBnO7Kbldz2ytugNUUddpPrnr9GFyV3U98k5yVDH/15OMX9Ay7xPFvDe2jf1sMX9Ry7xPF/N2aRv76faxxfx3FfO3b5u/BnBj+/hJ3uMdwN0ltsnHi3W+rW3+/yvmf3SQ+sid9eW6pd2em2R771vM+0SH5czq8doy9S4oxv4WWLFl/uPJgfTVDmPnd6jpTqbQvQT+j9zlX614vHWxruPJ/0sxp5i/QzF/py7vv1udfb/fPmp+HLnrf3OxzP37eM3EPtDvNG/QbVqsr9v/DF4M/HmS575RrH/FQcaPctkzZcLc7jTe3E7mdrf320fN5naX8U5TmzC3O403t5O53e399lFz7XKbmmbrjLzrSouDgNuAW4GLgLeRO3L/2sdrXwX8mdxBbHVwh7E7kHfGr7bNP3SSZb+BfO7ibyJi3YmJfCjZj4CtinMCu1mD/AXt16vJ26K9I/154Pbi+WHWV9YuwF+Aj7Q/kXp3HAep93MppQdb1vEn8uFlG/VZ713AphGxWZ/j250JrABsVTzehryfHg48iqWd0ReRO6iLBlzPhKm+X8id+xWAiXPzNu/jNbcALy4x3VKinrJWJf+y6WRJy5hBxo9y2TONub2Uub0sc9vcbmW21oe5vZS5vSxzuxm5XctsnemnrhxDPkdwBfJO8j7gCfR3cZkNgfNSSg+3zkwp3RwRi9vGbkDuWj/SNvbWDmMhH0q2CjkIJ7Mu+Zy+ydxN/oL2awPg/JTSQ201PhQRV7D08Khh1VfWRsCFKaVBLvwzSL3XdhhzB/l8vn7sTe4GXxwR15IPmzwVOLWPXxSQgxdy4P6w+HkW+RDHvxSPf1H8/F1Kqeuhn32Y0vuNiFeRfznvm1I6NyJuBXr+0ik+zx+XKXSE7gcePclzK7eMGWT8KJc905jby9Zobi9lbpvbrczW+jC3l63R3F7K3G5GbtcyW2d6o+OqlNLEh356RJwDnEM+/+p1fbw+TTI/plhXkA8BeneXMd1CBPKhnltHxIapuHL5EA2jvkFMtr17GaTehzuO6vOzTSl9J/J9tF9OPpd0O+DNwNkRsV1rN3eS198SEZeRz9dclXyv7HeklB6JiJ8C20bE0cA/km/fNFUDv9/ivMUjyRdFmqjlImB+RKzY7b1GxGzyoaz9uq39j50hugl4ekSslFJq70o/Hri97b2UGT/KZc805vZgzO0ezG1zu8N4DYe5PRhzuwdzuza5XctsnemNjmUUnbHjgd0i4oiU0rldhl8LPDUiZrfuCBGxPvnCO62uB54SEbNau4sR8WhgbodlX0XeGc/ssxvZybfJ55m9hXz1416uBZ4WEXNau8wRMQd4Kst2IIdRX1lXAptM8gXqZVT1dv1FUHR9TwBOiIgAPkU+93FH8v9s9HIm+WJDO5AP+/tJMf8n5MMwX0YOxjM7vrrPOofgk+QLf72y5btwEfmXzcbFvyfzRHpcQKnNxJXTR+E84CXAc8gXMwMgIlYGnkG+GNqg40e57BnN3Da3SzK3M3PbbK2MuW1ul2RuZ3XO7Vpm60y/RkcnHyV325Y7N63Nd8g7225t89/XYeyp5Cvxvr5t/j6TLPs44LFM0hGNiMf0qA3yxXWuAPaJiB0nWc6zIuJtxcNTyOH0lrZhexTzTx5yfWWdCKwFfLDD+np1QUdV773Fz2WuYhwRsyNfJfzvUkoTFyBabnwXZ5K/owcAN6aUrmmZvxKwH/AQLYFSps5hiIjnAXsCh6aULmx5aiJse503WKdzvb9O/iW1d9v8PcjnFZ44hfGjXLbMbXO7f+a2uW221oO5vSxze3Lmdv1zu5bZ6hEdbVJKV0fE14BdIuIFKaXJdupPA/8BHBsRzyLfL3o+8HzyxYRaHVyM/VJEPId8oZ6tyLcvup3lO4CHk3eyQyJiG/IX7W7gH8i3OlpCvihOt/dxf0S8Evg+cEpEnEG+ENAd5CB9EbB98T4m3s+/AUdFxBbkkHgm+fCvK1rGDaW+ARxO7rR+MCK2JN/yagmwKfA0cjez22tHUe955AsT7R8Ra5HvE30deXvdHBHfJW/HW8md0b3I5/ud2ufyzyqWvwn5PugApJQujYhbgKcDv0jFbcjK1plS+lWfdXQU+TZix5JvnXVQ29N9Be90nOsdEbuy9NzH9YAVI2LiF/gNKaXji1oujoijgLdHxEnAaeRt/07gp8BX2mrve/woly1zG3O7DHPb3DZba8DcNrdLMLdrntu1zdY0zbd5qcPE0lsC7TPJ85uQu8xntY1f0DbuH4Bvkb/E95C/UE8mH+qzqG3sBuQrTN9TjP9OMe92OtwWidyEeif5S3NfMV1F7oi9pMR7XZV8u6dzyF/6v5GvXv19YFdgdsvY9YDPAX8sxv2RfAuwdQetjyHd7qqYvzKwP/mX3BJgcbH+t/Xx2n7rXcDkt4ZaRNst7oDdyfd0f7B43ULyYW+fJN9L+g7yVYivJ98vfaOS++pviuXu2jb/xGL+x/rcdsvVOcj7bXt+f3Kgv7DDcysV+9D3a/B9X8Tkt89q/57OBt5D/uX5APAn8nmQq0+y7L7Hj3LZM2HC3Da3ze2Fg7zftufNbbN1Oj/HiX3b3Da3ze3m53btsjWKwlSBiFiHHLyfTyntWXU9kqTuzG1JGi/mtjQzeY2OaRKd73M9cX7hj6azFklSb+a2JI0Xc1vSBI/omCYRsQi4ATiffGjPtsArgXOBrdPobr0mSRqAuS1J48XcljTBRsc0iYj3kK8YPQ9YhXw+3knAQan3xW0kSdPM3Jak8WJuS5pgo0OSJEmSJDWG1+gAIuKGiLih6jokaSYziyWpWuawpKaYU3UBNbHmmmuuuSbL319bkqoWVRcwjcxiSXU1U7LYHJZUV6Vy2CM6JEmSJElSY9jokCRJkiRJjWGjQ5IkSZIkNYaNDkmSJEmS1Bg2OiRJkiRJUmPY6JAkSZIkSY1ho0OSJEmSJDWGjQ5JkiRJktQYNjokSZIkSVJjzKm6AEmS1DzHfv+CoS1rj1dsMbRlSZKk5vOIDkmSJEmS1Bg2OiRJkiRJUmPY6JAkSZIkSY1ho0OSJEmSJDWGjQ5JkiRJktQYNjokSZIkSVJj2OiQJEmSJEmNYaNDkiRJkiQ1xpyqC5AkSZIkDd+x379gKMvZ4xVbDGU5Vdvmk98eynLO3G/noSxHo+MRHZIkSZIkqTFsdEiSJEmSpMaw0SFJkiRJkhrDRockSZIkSWoMGx2SJEmSJKkxbHRIkiRJkqTGsNEhSZIkSZIaw0aHJEmSJElqDBsdkiRJkiSpMWx0SJIkSZKkxrDRIUmSJEmSGsNGhyRJkiRJagwbHZIkSZIkqTFsdEiSJEmSpMYYu0ZHRGwbEQsj4oqIuD8i/hgRJ0XE5lXXJkmSJEmSqjV2jQ5gT+AfgM8ALwPeXTw+LyKeV2VhkiRJkiSpWnOqLmAA/y+ldGvrjIg4A7gOeC+wcyVVSZIkSZKkyo3dER3tTY5i3mLgKuAJ01+RJEmSJEmqi7FrdHQSEesBmwGXVF2LJEmSJEmqzjieurKMiAjgGHLT5tBJxizusZg1h12XJGlZZrEkVcscljRTjH2jAzgEeBXwxpTSZVUXI0mSJEmSqjPWjY6I+DjwHuC/UkoLJxuXUprbYzmLsYMtSSNlFktStcxhSTPF2F6jIyI+AnwA2DeldETV9UiSJEmSpOqNZaMjIg4APgR8KKV0SNX1SJIkSZKkehi7U1ci4j3AgcD3gB9HxPNann4gpfTbSgqTJEmSJEmVG7tGB7BD8fOVxdTqBmDetFYjSZIkSZJqY+waHSml+VXXIEmSJEmS6mksr9EhSZIkSZLUiY0OSZIkSZLUGDY6JEmSJElSY9jokCRJkiRJjWGjQ5IkSZIkNYaNDkmSJEmS1Bg2OiRJkiRJUmPY6JAkSZIkSY1ho0OSJEmSJDWGjQ5JkiRJktQYNjokSZIkSVJj2OiQJEmSJEmNYaNDkiRJkiQ1ho0OSZIkSZLUGDY6JEmSJElSY9jokCRJkiRJjWGjQ5IkSZIkNcacqguYSW69efFQlvPo9ecOZTmSJEmSJDWNR3RIkiRJkqTGsNEhSZIkSZIao1SjIyJ2i4h5XZ6fFxG7TbUoSZIkSZKkQZS9RseXgF2B6yd5/rnFmOOmUJMkSZKkhvK6dZJGreypK9Hj+RWARwasRZIkSZIkaUoGuUZH6jQzIuYCrwBunlJFkiRJkiRJA+rZ6IiIAyLi4Yh4mNzkOGHicesE3AG8FvjaiGuWJEmSJEnqqJ9rdFxIvuZGALsBZwPXto1JwL3AL4GvDrNASZIkSZKkfvVsdKSUvgN8ByAingR8LKX0k1EXJkmSJEmSVFapu66klF40qkIkSZIkSZKmquztZQGIiFWBecA6dLgTS0rpZ1MrS5IkSZIkjcpjXn/QUJbz568eMJTlDFOpRkdErAYcBrxxktcG+Xods6demiRJkiRJUjllj+j4LPBm4DTgTPKdViRJkiRJkmqhbKPjVcBXU0q7jKKYOrhv8f1DW9Zqc1cd2rIkSZIkSVJvs0qOXwVYNII6JEmSJEmSpqzsER3nAxuNohBJkiTNHE2+CJ6a6VcXXD+U5Tx3i3lDWY6kyZVtdLwfODUivplSOm8UBUnDcOz3LxjKcvZ4xRZDWY4kSZIkaXqUbXS8Ffgj8IuI+AVwLfBw25iUUnrzMIqTJEn+L6IkSVIZZRsdC1r+/S/F1C6R78wiSZIkSZI0rUo1OlJKZS9eKqlPC47+4dCWtXDP7Ye2LEmSJEkaJzYuJEmSJElSY5Q9dQWAiFgNeD7wGODHKaU/D7UqSZIkSZKkAZQ+oiMi9gL+BJwBHAdsWsxfLyKWRMRbh1uiJEmSJElSf0o1OiJiZ+Ao4CzgLUBMPJdSug34AbDjMAuUJEmSJEnqV9lTV94LnJVSenVErAP8X9vz5wN7DKUySZJG4NabFw9lOY9ef+5QliNJkqThKnvqyubAyV2evxl49ODlSJIkSZIkDa5so+PhHq95HHDf4OVIkiRJkiQNrmyj43fA9p2eiIhZwL8B5021KEmSJEmSpEGUbXQcCbwsIj4KrD2xjIh4GvBN8h1YjhhifZIkSZIkSX0rdTHSlNLXI2JzYH9gv2L2D8h3XwnggJTS6cMtUZIkSZIkqT9l77pCSumDEXESsAuwMbnBcRVwfErp/CHX11FErA58gnyqzFzg98BHUkrfnY71S5IkSZKkeird6ABIKV0AXDDkWso4GdgC2Be4DlgAnBwRO6SUTquwLklj4DGvP2goy/nzVw8YynIkSZIkDc9AjY4qRcTLge2AnVJKJxfzzgI2BA4DbHRIkiRJM9DVV90ylOU8ZaPHDmU5kqrRtdERER8GEvDxlNIjxeNeUkrpo0OprrNXA3cB32ldYUR8GTgmIp6eUrp0hOuXJEmSJEk11euIjgPJjY6DgQeLx70kYJSNjs2AS1NKj7TNv6j1+RGuX5IkSZIk1VSvRscGACmlB1sfV2wd4MoO8+9seX4ZEbG4xzLXnGpRkqTuzGLNdNt88ttDWc6Z++08lOVo5jGHNUwHfO2coSznoNdtNZTlSK0ipVR1DaVExJXAFSmlHdrmb0RugOyVUjq67bmeob7mmmuyeHGvYfU1yvMRf3XB9UNZ9nO3mLfM41N/dvlQlrvD1hsPZTn9GsdQH+Uf1xu/98tDWfblh+y+3LzH7XHEUJZ907HvHMpy+rH2y/YaynLuPP1/J/4ZQ1lgDZTJ4vsW3z+Uda42d9WhLKcf43pe+Lhm8agsOPqHQ1nOwj23H8pyqjaOOQxm8WSqyGGY3iweR+bwsoaVwzC9WezfxMuqOodLXYw0IuYAq6aU7p7k+TWA+1NKD5VZbkl30OGoDWDt4ued7U+klOZ2W2AR+nawJWmEmp7FXrhOUt01PYclacKskuMPA87v8vx55Ot5jNLvgU0ior32zYufl4x4/ZIkSZIkqabKNjq2B7odA/9t4GWDl9OXk4G5wA5t83cjn9LihUglSZIkSZqhSp26AjwRuKbL89cWY0bpNOAs4AsRsQ5wHbA7sBWw44jXLUmSJEmSaqxso+NBYP0uzz8WaL/t61CllFJEvAr4RDHNJd9OdqeU0qmjXLckSZIkSaq3sqeu/BZ4bUSs2P5EMe/fgYuGUVg3KaW7U0pvTyk9NqW0ckppi5TSKaNeryRJkiRJqreyjY6jgE2B70fEsyNixWJ6NvA94OnAkcMuUpIkSZIkqR+lTl1JKX07Ij4J7Af8CkjFNIt8X9uDU0pfH3qVaqSm3OtbkiRJklQfZa/RQUpp/4g4BXgD8BRyg+MK4CsppfOGXJ8kSZIkSVLfSjc6AIqGhk0NSZIkSZJUKwM1OiRJGqXV5q5adQmSJEkaU10bHRHxYfI1OD6eUnqkeNxLSil9dCjVSTV00Ou2qroESZIkSdIkeh3RcSC50XEw8GDxuJcE2OiQJEmSJEnTrlejYwOAlNKDrY8lSZIkSZLqqFejY3fgpJbHCbgtpfTX0ZUkSZIkSZI0mH5OXbkauKR4fB2wK/CVEdakmnnuFvOqLkGSJEmSpL70anQsBua2PI4R1iJJy7jp2HdWXYIkSZKkMdOr0fFbYN+IWAH4SzHvBRHR9XUppeOGUZwkSZIkSVIZvRod7yZfo+MzxeME/GcxTSYBNjokSZIkSdK063Vkxu8i4qnAhsD6wCLg48CPR1+aJEmSJElSOV0bHRGxNXBZSukq4KqI+CmwKKX002mpTpIkSZIkqYRZPZ4/C3hxy+N5wGojq0aSJEmSJGkKejU6HgBWann8JGD10ZUjSZIkSZI0uF4XI70S2D0iLmDpXVfWiYh/6PailNKNwyhOkiRJkiSpjF6Njo8BXwEuKB4n4LPF1M3sKdalkp6y0WOrLkGSJEmSpMr1uuvKtyLid8B88l1XDgBOAS4afWmSJEmSJEnl9Dqig4k7rgBExIHAt1NKXxlxXZIkSZIkSaX1bHS0Sin1unipJEmSpDG32txVqy5BkgY2UOMiIraOiI9FxLERsXExb/Vi/tzhlihJkiRJktSfUkd0RMRs8sVJXwME+eKkXwUuBx4iX7/jUOATwy1T0lScud/OVZcwY9x5+v9WXYIkSZI0o5U9ouN9wM7Au4FNyM0OAFJKS4CTgZcPrTpJkiRJkqQSyjY6dgOOSykdDtze4fnLgCdPuSpJkiRJkqQBlDp1BZgHHNbl+cXAWgNXI0mSptUOW29cdQmSJElDVfaIjnuAtbs8/xTgtsHLkSRJkiRJGlzZRsc5wBsiItqfiIi1gDcBZw2jMEmSJEmSpLLKNjo+DmwEnAm8spj3TxHxn8AFwGrAp4ZXniRJkiRJUv9KXaMjpXR+ROwEfAH4UjH7UPLdV24FXp1SunS4JUqqs8sP2b3qEiRJkiTp78pejJSU0mkRMQ94CbAxuclxFfDDlNL9Q61OkiRJkiSphNKNDoCU0gPAqcUkSZIkSZJUCwM1OiJiDWA7YMNi1rXAj1JK9wyrMEmSJEmSpLJKNzoi4i3AYcDq5NNWABJwb0S8O6X0hSHWJ0mSJEmS1LdSjY6I+FfgGPIRHB8GLime2hR4B3BMRNyaUvKUFkmSJEmSNO3KHtGxL3AZ8NyU0r0t838SEV8Cfgm8D6/dIUmSJEmSKjCr5Ph/Aha2NTkAKK7P8eVijCRJkiRJ0rQr2+iApdfl6CQNWogkSZIkSdJUlW10/A7YPSJWa38iIlYHFhRjJEmSJEmSpl3Za3QcCpwEXBARRwCXFvMnLkb6FGCn4ZUnSZIkSZLUv1KNjpTSKRHxduBg4H9YeqpKAPcBb08pfWe4JUqSJEmSJPWn7BEdpJQ+FxFfAV4MbEBuclwD/CildNeQ65MkSZIkSepb6UYHQEppMfDNIdciSZIkSZI0JT0vRhoRsyPiUxGxZ49xe0XEJyKi211ZJEmSJEmSRqafu668AXgvcF6Pcb8G3ge8fqpFSZIkSZIkDaKfRsdrgR+nlH7TbVDx/A+x0SFJkiRJkirSzzU6ngUc1ufyzgLePXg5kiRJ1Vm45/ZVlyBJkqaon0bH2sCtfS7vtmK8JEmSxtxNx76z6hIkSSqtn1NX7gHW7XN56wD3Dl5ObxGxc0R8PSKujYi/RsR1EfHliJg3yvVKkiRJkqT666fR8XvgJX0u78XF+FHaF1gZ+AjwUuBA4J+BCyJigxGvW5IkSZIk1Vg/jY6TgO0iYsdugyLiX8mNjm8Po7Audkgp7ZhSWphS+mlK6cvkRsxc4O0jXrckSZIkSaqxfhodnweuBr4RER9vP0UkIuZFxMeAbwBXFuNHJqW03PVCUkrXAbcDTxjluiVJkiRJUr31vBhpSumvEfEK4HvAfsD7I+Ie4G7gUcAaQABXAK9MKS0ZYb0dRcRmwHrAJZM8v7jHItYcelGSpGWYxZJULXNY0kzRzxEdpJSuBp4B/BdwDvAQ8FjgYeDsYv4WKaVrRlTnpCJiJeALwB3A0dO9fkmSJEmSVB/93F4WgOJIjf8ppqGIiPnAWX0OXy+ldHvb62cDx5GbMK9MKd3W6YUppbk96liMHWxJGimzWJKqZQ5Lmin6bnSMyOXAG/sce0/rg4iYBXwJ2An495TSj4ZcmyRJkiSN1A5bb1x1Caq5m459Z9UljJ1KGx0ppVuAhWVfVzQ5vgj8B/CGlNJJQy5NkiRJkqS+XH7I7lWXoBZVH9FRWkQEcCywK/DGlNLXKi5JkiRJkiTVxNg1OoAjgDeRmx1XRsTzWp67O6V0aTVlSZIkSZKkqo1jo2OH4ucexdTqp8D8aa1GkiRJkiTVxtg1OlJK86quQZIkSZIk1dOsqguQJEmSJEkaFhsdkiRJkiSpMWx0SJIkSZKkxrDRIUmSJEmSGsNGhyRJkiRJagwbHZIkSZIkqTFsdEiSJEmSpMaw0SFJkiRJkhrDRockSZIkSWoMGx2SJEmSJKkxbHRIkiRJkqTGsNEhSZIkSZIaw0aHJEmSJElqDBsdkiRJkiSpMWx0SJIkSZKkxrDRIUmSJEmSGsNGhyRJkiRJagwbHZIkSZIkqTFsdEiSJEmSpMaw0SFJkiRJkhrDRockSZIkSWoMGx2SJEmSJKkxbHRIkiRJkqTGsNEhSZIkSZIaw0aHJEmSJElqDBsdkiRJkiSpMeZUXYAkSZIkSQAL99y+6hLUAB7RIUmSJEmSGsNGhyRJkiRJagwbHZIkSZIkqTFsdEiSJEmSpJ6x3h0AACAASURBVMaw0SFJkiRJkhrDRockSZIkSWoMGx2SJEmSJKkxbHRIkiRJkqTGsNEhSZIkSZIaw0aHJEmSJElqDBsdkiRJkiSpMWx0SJIkSZKkxrDRIUmSJEmSGsNGhyRJkiRJagwbHZIkSZIkqTFsdEiSJEmSpMaw0SFJkiRJkhrDRockSZIkSWqMOVUXIEmSJA3Tnaf/b9UlSJIqZKNDkiRJkiQNTdUNZ09dkSRJkiRJjTH2jY6IWBgRKSJOqboWSZIkSZJUrbFudETEi4HXAHdXXYskSZIkSare2DY6ImJ14FjgAOAvFZcjSZIkSZJqYGwbHcAngTuBz1ZdiCRJkiRJqoexvOtKRPwz8J/A81NKD0dE1SVJkiRJkqQaGLtGR0SsBHwBODKl9Js+X7O4x5A1p1yYJKkrs1iSqmUOS5opKj11JSLmF3dM6Wdat3jZAcAqwIcqLF2SJEmSJNVQpJSqW3nEY4GX9jn8q8CGwEXA7sBpLc9dBFwM7ALcn1J6sGQdjwCx5po2sSXVy1133XVjSulJVdcxHcxiSXU1U7LYHJZUV2VzuNJGR1kR8Srg5B7D9kopHV1yuQ+Rj27p5za1E8l/V5l11MA41j2ONcN41j2ONcN41l225rtmwh/XMCOyeBxrhvGsexxrhvGsexxrBrO4oxmQwzCedY9jzTCedY9jzTCedY80h8et0bEusFmHp74GXA18ELgypXTTCGtYDJBSmjuqdYzCONY9jjXDeNY9jjXDeNY9jjXX0Thux3GsGcaz7nGsGcaz7nGsGca37joZ1204jnWPY80wnnWPY80wnnWPuuaxuhhpSul2YFH7/IhYAtyeUlruOUmSJEmSNHNUejFSSZIkSZKkYRqrIzomk1KaV3UNkiRJkiSpeh7RIUmSJEmSGsNGhyRJkiRJagwbHZIkSZIkqTFsdEiSJEmSpMaIlFLVNUiSJEmSJA2FR3RIkiRJkqTGsNEhSZIkSZIaw0aHJEmSJElqDBsdbSLiSxFxf0Q8tcNzn4mIhyJiy+LxxyPi9Ii4NSJSRBw47QUvra2vuiNi3Yj4VkRcHRH3RsTdEfHbiHh7RMyuY83F4zTJtGcda46IBV1qThHxujrWXTx+QkScGBF3RMRfI+L8iHh1jerr63sXETtExAkRcVlEPBwR149J3R8otvkdEfFARFwXEcdExBOHXX+djWMWm8P1q7tOWVz3HB6gxlpksTk8OuZw/eouHtcii8cxh8vUXTz2b+Jq655aFqeUnFomYE3gj8AvgNkt87cCHgY+3jLvXuDnwDFAAg6se93AE4ATgDcD2wEvBQ4v6j+6jjUX8xLwNeB5bdOj61gzsF6HWp8H/Ar4KzC3pnWvBdwA/AFYUOwfxwOPAK+pur5iXl/fO+ALwOXAicAVwPVjUvengfcDOwDzgb2APwE3AWtN535T5TSKbVuXmjGHp3Nb1yaLS9RcSQ4PsD/UIotHVLM5PKJtW5eaqVEOD7Cta5HFJbZ1bXK4ZN3+TVx93VPK4mnbqcZpAl5WbPR9i8erAVcDFwMrtoybVfyc2+1Dqlvdk7z2a8ADwJw61lyM+WzV+8ZUtjPwaOBB4MS61g18oAjwZ7a99izgxol9vuJ9oa/vXWutwCmjCPVR1D3JOl5avGa3KvadqqZxzGJzuH51d3hdZVlc9xwuuT/UJovN4fHZH+pU8ySvrSSHS27r2mTxOOZwv3VXmcXjmMOjqHuSdfSdxdO+Y43LBHwRWAI8HTgS+BuwxSRjKw/1Qepue92R5K7qSP+AGrTmOoX6oNsZ2Kd4H9vWtW7gVOCPHV73rqL259dlu5b53o0y1EdZd8trnl285vVV7TtVTeOYxeZwveru8JpKs7juOVx2u9Yli83h8dsf6lJz2+sqy+F+665bFo9jDvdTd9VZPI45PMq6W17TdxZXsmONw0Q+/OYPwDXkbt5BXcbWItTL1A0EMId8WNZryYcQfaSuNRfb987il88S8uFur637dm57ze+B64Coa93AD4FrOrzubcVnsEddtmvNQn3odRffz1WAZwDnAJcBq1W171Q1jWMWm8P1qrvDayrN4rrncNntWpcsNofHb3+oQ811yuF+665bFo9jDvdTd9VZPI45PKq6B83iSnascZmAtxQb/wpghWHsXHWpG3h7MSYVO+HH6lwz+TzK/wBeUPwiWlSM/6+61tw29nnF2A/Xef8APgM8BDy+bf6JxWv2q8t2rVOoD7tuYPWW72cCfgmsX/W+U9U0jllsDten7raxtcjiuudwme1apyw2h8drf6hDzXXL4X7qrmMWj2MO96q7Dlk8jjk87LqnksXedWUSETGH3LF7BHgSsNwVZOuoRN1fB7YEXgJ8CtgnIv5nWops00/NKaU3pJS+klI6O6X0DWAb4GzgYxGxyrQWzED7x5uKsQtHW1l3fdR9DDnUvxIRm0bEOhHxDuA1xfOPVFxfLY2g7vvJ389/IV8obS1gUUSsP8Xljp1x3CfM4ekzjllc9xzus8baMYdHp+H7Q21yGMYzi8cxh6H+WTyO3zuoWRZPd+dsXCbgw8UHtCP5gjO/puUKsmW7UXWsu+117ynewzNHXeMQa35rUfOWda4ZWBW4CzhjHPYP8oWEbmRp5/RGlv6Px65V19cytjbd61HnBfA48sXRDq96H5ruaRyz2ByuZ911yeK65/AA27UWWWwOj9/+UJea215XWQ5PsW7/Jh5B3VVm8Tjm8CjrbnlN31lc6Q5W1wn4R/KVgI8sHm9PyxVkh/Eh1aHutte+oBj772NU857F2GfVuWZgtyq27VTqJp+zuhGwCTCbfIjkI8AGdaiveL4WoT5deUE+1/H0Kveh6Z7GMYvN4frWXYcsrnsOD7hdK89ic3g894e61Nz22kpyeAh1+zfxiOquIovHMYdHXXfb6/rK4sp2sLpOwArAb4FrabnICfn+w38FNh7Wh1R13W2v//B0B+RUagZmkQ/TuxtYuc41k8+dvANYaRz3D2Alcjf2e3Wqrw6hPl15AWxIPnzyf6rah6Z7GscsNofrva2rzuK65/AUtmulWWwOj+/+UJea214/7Tk81bqryuJxzOEhbGv/Jq6o7pbX9J3FlexgdZ6AA8hduhe1zZ+4guy5LL337wvJ52lNdCe/UTx+DbBqHesm387pC8AuwHzyYUVHFTvMt2pc87HA64uaXwf8tNjmb6tjzS3zNyzGV/qHUYltPQs4Atip2NZvBi4iH3r2xKrrK+b19b0jnxc4Mf/XwK0tj59ex7qLms8mHxb5UmA7YO9i+98GzKtyP6rjPltmn6hLzZjD07p/FPMrz+IS27qSHC67Xfv93jHiLB52zZjDI90f6lIzNcrhAequRRaX2T+K+ZXncMlt7d/EFdbNELK4sp2sjhPwT+TDbY6a5PmXFR/Ge4rHi1j2KrCtU8+NX0XdxU7yA+Dm4jX3AOcB/wXMqWnNOxQ7+m3kezEvBn4C7FDn/aOY91EqPNdzgG09C/gucFPxmj8CRwOPrUN9xeO+vnfAgi7jDqxj3eRfBguBK8m3uHuA3Bn/PPCkqvahOu+zZfaJutSMOTyt+0cxr9IsLrmtpz2HB9mu/X7vGGEWj6JmzOGR7g91qZma5PAAddcii8vuH8U8/yYe8nbt93vHDP2bOIoFSZIkSZIkjT1vLytJkiRJkhrDRockSZIkSWoMGx2SJEmSJKkxbHRIkiRJkqTGsNEhSZIkSZIaw0aHJEmSJElqDBsdkiRJkiSpMWx0aCxExPyISG3TvRFxQUS8KyLmVF2jJDWdWSxJ1TKHpf74RdC4+SpwGhDAY4HdgP8GNgHeWmFdkjSTmMWSVC1zWOoiUkpV1yD1FBHzgbOA96aUDm2ZvxpwOfB44DEppdsqqG02sFJK6f7pXvewRMSjUkr3VF2HpHozi0fLLJbUizk8WuZwc3jqisZaSuk+4JfkbvaTW5+LiGdHxMkRcXtEPBARV0TE/p0O6YuInSPidxGxJCJujIgDImK74nDABS3jFhTztouID0XENcAS4LVl1xsRm0bENyPiT8W4WyLirIh4RcuYlSPiwGIZ90fE4oi4OCIO6fAe3lIctvjXiLgrIs6IiK06jEsRsTAito2IcyLiXuDUEptdkpZhFi+zPLNY0rQzh5dZnjksT11RI0yE+Z0TMyLi5cDJwNXAYcVzzwc+AjwD+LeWsf9OPvzvGuAg4CFgd2CHLus8FFgBOBa4G7iizHojYh3gzGJZRwM3AOsCzwaeC3y/eO4o4E3AccBngNnARsA2rcVExMHAvsCvgQ8AjyIftnhWROyYUjqtrf5nAzsX9X+5y/uUpH6ZxWaxpGqZw+awJqSUnJxqPwHzgQR8mBx+6wGbk0MvAb9uGbsycAvwM2BO23LeVYyfXzyeA/wJ+DOwVsu41YFri7ELWuYvKOZdAazatuwy6/3X4vFre7zvO4HTeox5GvAIcA6wYsv8xwGLgeuB2S3zUzFtV/Xn6uTkNF6TWdx1jFns5OQ08skc7jrGHHb6++SpKxo3BwG3AbcCFwFvA04ih+SEFwOPAb4EzI2IdScm8kWbAF5S/HwWOfwWppT+MrGAlNK95K7yZP43LX/+YZn13lX8fFlErNFlPXcBm0bEZl3G7Eg+TPHTKaUHW97DTcBC4EnAM9te87uU0o+7LFOSujGLl2cWS5pO5vDyzGH9naeuaNwcA3yTfIjc5sD7gCeQzwmcsEnx84tdlvOY4ucGxc8rOozpNG/ClR3m9b3elNJPI+I4cjd8l4g4D/gx8PWU0qUt4/cGjgcujohryRefOhU4NaX0SNt7+H2H9V1S/NwQOL9H/ZLUL7PYLJZULXPYHFYXNjo0bq5q6bqeHhHnkA9POxp4XTE/ip/vBS6cZDk3tY0tq9PVpMusl5TS7sUFlF4ObAW8B9g/IvZOKR1ZjPlORMwrxrwQ2A54M3B2RGxXdKsHeQ9jezVsSbVgFpvFkqplDpvD6sJGh8ZaSunciDge2C0ijkgpnQtcVTx9Xx+Hol1X/Hxah+c6zeumzHoBSCldQu4wfzoi5gK/Aj4VEUelVJw8mNKdwAnACRERwKfIF1nakdzJv6ZY3KYt/57w9OLntSXfiyT1zSw2iyVVyxw2h7Usr9GhJvgo8DD5Ks4APySfr/j+iFi7fXBErBIRjyoeng/cDCyIiLVaxqwO7Fmyjr7XGxFrR8Qy37+U0mLyL5lVgZUjYnYR9K1jEvDb4uHEOr5LvpDSeyNihZb1rQ+8kXz16t8iSaNlFpvFkqplDpvDKnhEh8ZeSunqiPga+by+F6SUzo6I3YBTgCsi4ovkW1vNBTYGdgJeDSxKKT0UEfsAJwK/jogvkG+ltQC4g3yuX+qzjvv6XS+wG/CuiJi47dbfyIfhbQ98I6X01yLQb46I75JD+dainr2Av1Dc5zuldEVxuN++wM8i4ussvZXW6sAuKaWHy25XSSrDLDaLJVXLHDaH1aLq2744OfUzsfRWWvtM8vwm5A72WS3zNiMf3vYn4EHy7bLOBT4ErN32+teSr1j9AHAjcAA5gJe53RVLb6U1v0utPddLvn/4l8mBfh/5vuO/I5+TuFIxZkXgk+T7gN9R1HY9+cJOG3VY7x7k8F9SLO9HwAs6jEvkK2pX/rk6OTmN12QWm8VOTk7VTuawOezU3xTFhyypTUS8BzgUeH5K6ZdV1yNJM5FZLEnVMoc1jmx0aMaLiBWBh1PLoWzF+YgXAWsAj0st9+KWJA2fWSxJ1TKH1SReo0PK99M+vTin8TpgfWB3inP/DHRJmhZmsSRVyxxWY9jokOA24JfALsCjyRdeuhh4f0rpG1UWJkkziFksSdUyh9UYnroiSZIkSZIaY1bvIZIkSZIkSePBRockSZIkSWoMGx2SJEmSJKkxbHRIkiRJkqTGsNEhSZIkSZIaw0aHJEmSJElqDBsdkiRJkiSpMWx0SJIkSZKkxrDRIUmSJEmSGsNGhyRJkiRJagwbHZIkSZIkqTFsdEiSJEmSpMaw0SFJkiRJkhrDRockSZIkSWoMGx2SJEmSJKkxbHRIkiRJkqTGsNEhSZIkSZIaw0aHJEmSJElqDBsdkiRJkiSpMWx0SJIkSZKkxrDRoZGLiA0i4pSIuC0iUkQsnGx+RMwv/r1ggPUM/NqmKLsNImJBMX7+ENZ9bkQ8EBG/jIh5U12epOqY29PH3JY0DOb29DG3x8OMbHS07Jz7VF3LKEXEqhGxd0ScHRF3RsTfIuLPEXFa8YWbM02lLAReCBwM7Ap8vsf8sRARz4iIA+seMtNY52HAccBzgVp8tyJiv4j4ZkRcW3znr+8xflZEvCsiLo+IJRHxh4g4LCJWm8rYUS97JjC3ze1hMLeXY27P8GwdJXPb3B4Gc3s5dcztemZrSmnGTcB8IAH7VF3LCN/jU4Arivf5I+C9wBuB9xWPE/DpaahjJeAR4Ig+588CVgZmD7CugV874HtbUGzH+VV/3t22Qbc6h/0egDnAPcAvqt4WRT0JuKPY5+8Eru8x/vDiNScBewD/DfwNOBOYNejYUS97Jkzmtrk9pPdmbi+/PHN7BmfriD9Lc9vcHsZ7M7eXX17dcruW2TpdHUZNo4hYBfgesCGwc0rppLYhB0fElsCW01DOY4Ag/7HSc35K6RFgySArmsprm6LqbZBSeigiLgE2i4hIRfpV6MkppWsBirpWn2xgRGwKvAM4KaW0c8v864AjgNcBXyk7dtTLVjOY2zNX1dvA3DZbNRhze+aqehvUKbdrna1Vd4Aq6jrNp48OM/Ao4GPAr4DbgQeAq4FPAau2jV0ZOJDc1b0fWAxcDBxSZkwxbl3gKOAPwIPFz6OAdfp8f+8o3t+nSm6XvtdL7hB/APg9+Yu+GDgVeGbLmIVFHe3T9ZPMn9/y2SxoW9+KwL7AhcW2uws4H3h7h8+1/bU9ay3GLShevw35ULBris/8SmD3trEHTvIeFpb5rDts1ycVyzmwbf4Zxfy92+b/Crh0sm3QR519v+c+96EAflMsc4Oqv+tttV1Cl/8ZJH/XE/CCDt/t+4DTBhk76mXPlAlze7LXmdvmtrltttZywtye7HXmtrndmNymxtnqER3dPR54C/BtcifqIfI5bvsCzwS2bxl7FPAm8jlTnwFmAxuRd+i+x0TEmsC55EPhvghcUKxrL2CbiHhOSumeHnW/pvh5TL9vtMx6I2IF4AfAPwPHA0cCa5IPVfp5RGydUjqffA7ghcV7PZl8OBPkQH9mh/mXAZt0qG1F4IfkUDkDOIEcoJsDOxXrn+x99Vtrq08AqxT1P1Bsg4URcXVK6efFmJOA9YG3FuMvK+ZfU/zsZ39YTkrphqIDui05NCfe/7+QDz3cFvhsMX8N4FnA0V0W2avOMu+5H3sBWxT/3hy4brKBETELWLvEsu9MuYM+KluSt/GvW2emlJZExIUs+z8yZcaOetlalrltbpvb5rbZOl7MbXPb3B7f3K5vtlbZAaqw8zSf/jrMKwIrdJj/0eL1z2mZdyc9OlZ9jvl4sey3tc3/f8X8j/bx/u4A7i65TfpeL/CuYt72bWPXAG4EFrXMm0fnjulk8yc+mwUt8/Yt5n2iQ92zery2TK0LirG/BVZsmf94chh9tW0ZE+PnD/JZd/ks/o/c4V+teLx1sZ7jgbuBOcX8HYr5O/XYBt3qLPWee9T9OHLn/+Zimfv3GD+xD/Q7zRtke7asr9f/DF4M/HmS575R1LBi2bGjXvZMmTC3p7RezO328fMH+ay7fBbmdsNy22nqE+b2lNaLud0+fv4gn3WXz8LcHkJuU+NsnZF3XelXSunBlNLfACJiTkSsFRHrAj8uhjy3ZfhdwKYRsVmXRfYz5tXAbSzfHf48+XC+V/dR+hrkL2gZZdb7BuBy4DcRse7ERP5F9SNgq+K8xWHZBfgL8JH2J1LvjuMgtX4upfRgyzr+RD60bKMSNffzWU/mTGAFYKvi8TbAreQL/TyKpZ3RF5E7qIsGWEe7YbznI8l1T5yft3mP8bcALy4x3VKilkGsSv5l08mSljFlx4562Wphbk+6XnO7N3Pb3DZbK2BuT7pec7s3c7v63K5ttnrqSg8R8TZgT2BTlr8d71ot/96b3AG8OCKuBc4in5d2aks49DNmA+D8lNJDrStK+aIzV7D0MKVu7iZ/Qcsos95NyIdd3dZleeuSD5kbho2AC1NKg1z0Z5Bar+0w5g7y+Xz96ueznsyZxc9tyIcQblO8/gLyL6BtgF8UP3+XUmq/8NQgpvSeI+JV5F/O+6aUzo2IW4Guv3SKz/PH3cZMs/uBR0/y3MotY8qOHfWy1cbcNrcL5nYX5nbXsZpm5ra5XTC3u6hpbtc2W210dBER7ybfq/gM8lVjbyIf4vR48oV//h7EKaXvRL5v8svJ5xVuB7wZODsitiu61T3HDKn0S4CtI2LDVFy5fMiCfJjSu7uM6RZ0g0gDvm6QWh/usqy+TOWzTindEhGXkc/VXJX8PxnvSCk9EhE/BbaNiKOBfyTfvmkYBn7PxbmLR5IvijRRz0XA/IhYcbL3GhGzgfVK1HhbSmmyOofhJuDpEbFSSqm9M/144PaW91Jm7KiXrRbm9qTM7R7MbXPbbK2GuT0pc7sHc7sWuV3bbLXR0d2u5CsWv6y1KxgRL+00uOj0nQCcEBFBvlr0vsCOwDf7HHMt8LSImNPa7Y2IOcBT6dwJbPdt8nlmbyFf/bgfZdZ7FfkLc2Yf3dJhuBLYZJIvUC+jrLXrL4N+9ocuziRfaGgH8mF/Pynm/wQ4FHgZORTP7PjqEnUOwSfJty97ZUswXkT+ZbNx8e9OnkiXiyd1sAH5+zgq5wEvAZ4DnD0xMyJWBp4B/GzAsaNetpZlbpvb3Zjbmbndfayml7ltbndjbmd1ze3aZqvX6OjuYfJO+/cuWxFE728dFBGzI2Ju67yU0sQFZwDW7mdM8fMUclC8pa2WPYr5J/dR9/+Rb7O0T0Ts2GlARDyrOExwQpn1Hgc8lkm6thHxmD5qLONE8mGLH+ywrl4d0FHWem/xc5krGZf4rLs5k/z9PAC4MaV0Tcv8lYD9yFclP7vzy3vXOQwR8TzyoaaHppQubHlqImy7nTdYt3O9v07+vu/dNn8P8rmFJw44dtTL1rLMbXO7G3Pb3DZb68fcNre7Mbfrndu1zdaZfkTHtkW3qd3tKaWjgW+Ru2enR8RJ5IsO/Qfwt7bxjwJujojvkr9ct5K7YXuRz/E6tc8xAJ8G/g04KiK2KMY+k3wY1hXF812llO6PiFcC3wdOiYgzyBcCuoMcoi8i36qrdVll1ns4+YtwSERsQw6Du4F/IN+OaUmxjmE5nNxp/WBEbEk+tHEJ+TzOp5E7md1eO6pazyNfnGj/iFiLfK/o68jbq5/PupuzimVvQj5sE4CU0qURcQvwdOAXqfetzyatM6X0q37e5GQi30rsWPKtsw5qe7pn8E7Hud4RsStLz3tcD1gxIiZ+gd+QUjq+pZ6LI+Io4O3F9/008vZ/J/BT8i3vSo8d9bJnIHN7KXO7PHPb3DZbp5+5vZS5XZ65XePcrnW2pgpu9VL1xNJbAk02XV6Mm03u5F1NvprsDeQA2oSWWzWRD3X6JPn+wXcUY68n3x97o37HtNS3HvA54I/kkP8j+T7R65Z8n6uSb/d0DvkL/zfgz+RA3hWYPeh6yU2yd5K/1PcV01Xkrt1LWsbNa91Wfcyf+GwWtM1fGdgf+D05LBcX635bH6/tt9YFTH5bqEV0uMUdsDtwKflc0kQOyb4/6x6f32+KZe7aNv/EYv7Huuzb7dtguToHfc8tz+9PDvQXdnhupWIf+n7F3/VFTP49X9Rh/GzgPeRfng8AfyKfB7n6VMaOetkzYcLcNrfN7YWDvueW581ts3U6P8v5XT5Lc9vcnt/h/S7C3G5/fhxyu5bZGkVxkiRJkiRJY89rdEiSJEmS/n97dx5mW13eif77Ak5A5DgPbVq0NcEpAzGP2q02neAUJQ4kxrQGMIYEvXYGjWNawSlq1Gti9MaWa0TjnDiFKDEx4nQTp8Z2BlRQ2yEBJQcHRFr43T/2OqEoz6nau2rvWsP5fJ5nPcVew299a9Wq99TzsgaYDI0OAAAAYDI0OgAAAIDJ0OgAAAAAJkOjI0lVfamqvtR3DoD9mVoM0C91GJiKg/oOMBCHHXbYYYdl9tofgCGpvgPsILUYGKr9pRarw8BQLVSHXdEBAAAATIZGBwAAADAZGh0AAADAZGh0AAAAAJOh0QEAAABMhkYHAAAAMBkaHQAAAMBkaHQAAAAAk6HRAQAAAEzGQX0HAACm59S3n7W0sU6875FLGwsAmD5XdAAAAACTodEBAAAATIZGBwAAADAZGh0AAADAZGh0AAAAAJOh0QEAAABMhkYHAAAAMBkaHQAAAMBkHNR3AAAAAJbv1LeftZRxTrzvkUsZp28/9+w3LWWcdz/p2KWMw+q4ogMAAACYDI0OAAAAYDI0OgAAAIDJ0OgAAAAAJkOjAwAAAJgMjQ4AAABgMjQ6AAAAgMnQ6AAAAAAmQ6MDAAAAmAyNDgAAAGAyNDoAAACAydDoAAAAACZDowMAAACYDI0OAAAAYDJG1+ioqp+vqtOq6pyquqSqvlJVb66qO/SdDQAAAOjX6BodSU5K8u+TvDDJfZI8pvv8kaq6c5/BAAAAgH4d1HeALfi/WmsXrJ1RVX+X5Pwkj0tybC+pAAAAgN6N7oqO9U2Obt7uJJ9LcrOdTwQAAAAMxegaHXtTVTdIcvskn+o7CwAAANCfMd66chVVVUlellnT5vn7WGf3JsMctuxcAFyVWgzQL3UY2F+MvtGR5HlJHpDk4a21z/YdBgAAAOjPqBsdVfWsJI9N8juttdP2tV5rbdcm4+yODjbASqnFAP1Sh4H9xWif0VFVT0/y5CSPb629qO88AAAAQP9G2eioqpOTPCXJU1prz+s7DwAAADAMVNIK+QAAIABJREFUo7t1paoem+SUJH+T5F1Vdec1i7/fWvtYL8EAAACA3o2u0ZHkmO7r/bpprS8lOXxH0wAAAACDMbpGR2vtqL4zAAAAAMM0ymd0AAAAAOyNRgcAAAAwGRodAAAAwGRodAAAAACTodEBAAAATIZGBwAAADAZGh0AAADAZGh0AAAAAJOh0QEAAABMhkYHAAAAMBkaHQAAAMBkaHQAAAAAk6HRAQAAAEyGRgcAAAAwGRodAAAAwGRodAAAAACTodEBAAAATIZGBwAAADAZGh0AAADAZGh0AAAAAJOxUKOjqo6rqsM3WH54VR233VAAAAAAW3HQguu/IsmvJfniPpbfqVvnVdvIBAAATNQFX9+9lHFueJNdSxkHmJ5Fb12pTZZfLckVW8wCAAAAsC1beUZH29vMqtqV5L5Jvr6tRAAAAABbtOmtK1V1cpKndh9bkldX1as32OQFywgG+5sTXvrOpY112kn3WtpYAAAAYzLPMzr+V2bP3KgkxyV5f5Lz1q3TknwnyQeTvG6ZAQEAAADmtWmjo7X2tiRvS5KqunmSZ7bW/mHVwQAAAAAWtdBbV1pr/2VVQQAAAAC2a9HXyyZJqurgJIcnuV728iaW1tr7thcLAAAAWJUb/erTljLOv7zu5KWMs0wLNTqq6pDMHjb68H1sW5k9r+PA7UcDAAAAWMyiV3T8cZJHJHlHkncn+ebSEwEAAABs0aKNjgckeV1r7aGrCAMAAACwHQcsuP61krxnBTkAAAAAtm3RKzo+muTWqwgCsFOm/OAlgLFQixmbD531xaWMc6cjD1/KOMC+LdroeGKS06vqL1trH1lFIPYfp7/v7KWMc8zdj1jKOAAAAIzfoo2O30zylST/VFX/lOS8JJevW6e11h6xjHAAgP+LCACwiEUbHSes+e//1E3rtczezAIAAACwoxZqdLTWFn14KQAAAMCO0bgAAAAAJmPRW1eSJFV1SJK7JLlRkne11v5lqakAAAAAtmDhKzqq6pFJvprk75K8Ksntuvk3qKpLq+o3lxsRAAAAYD4LNTqq6tgkL0lyZpLfSFJ7lrXWLkzyt0nuv8yAAAAAAPNa9NaVxyU5s7X2wKq6XpL/d93yjyY5cSnJAGAFLvj67qWMc8Ob7FrKOAAALNeijY47JHnCBsu/nuSGW48zbf64ZoqOeNwrlzLO2c87finjAAAA+7dFn9Fx+Sbb3DTJd7ceBwAAAGDrFm10fDzJvfa2oKoOSPLLST6y3VAAAAAAW7HorSsvTvK6qnpGZm9cSZIDqurHk/xhZm9g2ejWFhi9k1//gaWM87SH3HUp4wAAAHClhRodrbU3VNUdkvxBkid1s/82s7evVJKTW2tnLDciAAAAwHwWvaIjrbX/XlVvTvLQJEdk1uD4XJK/aK19dMn5AAAAAOa2cKMjSVprZyU5a8lZ5lZVh2Z2q8wvJ9mV5NNJnt5a++u+MgEAAAD921KjYwDekuTIJI9Pcn6SE5K8paqOaa29YzsDf3f3JdtP1zlk18FLGwsAANjY5z/3z0sZ51a3vvFSxgH6sWGjo6qemqQleVZr7Yru82Zaa+0ZS0m390y/kOToJA9qrb2lm3dmklsmeUGSbTU6AAAAgPHa7IqOUzJrdDw3yWXd5820JCtrdCR5YJKLk7zt33bYWquqVyZ5WVXdtrX2mRXuHwAAABiozRodt0iS1tplaz/37PZJPtNau2Ld/E+sXb52QVXt3mTMw5aUDYB9UIvZ3/3cs9+0lHHe/aRjlzIO+x91mGU6+fUfWMo4T3vIXZcyDqxVrbW+Myykqs5Ncm5r7X7r5t86yblJHtVa+7N1yzYt6ocddlh2795steFa5f2IHzrri0sZ+05HHr6UceZx6tuX86zcE+975FLG6dtY/7i+6YkvWso4Xzv1t5cyzjyue59HLmWci874tzJWSxlwABapxct6XtJOPitprPeFn/6+s5cyzjF3P2Ip4/TthJe+cynjnHbSvZYyTt/GWIcTtXhf+qjDiefWbUYdvqpl1eFkZ2vxEY975VLGOft5x//QvDHW4r7r8EIPI62qg5Ic3Fr71j6WXzvJJa21Hywy7hZs1J35oWWttV0bDdYVfR1sgBWaei324Dpg6KZehwH2WPStKy9Icp8kP7aP5R9J8jdJHrudUJv4ZpLr7WX+dbuvF61w3zBKLnMGAAD2F4s2Ou6VZKNr4N+U5AFZbaPj00mOraoD1j2n4w7d10+tcN/7pZ285QQAAAC244AF1//RJF/YYPl53Tqr9JYku5Ics27+cUnO8cYVAAAA2H8tekXHZUlussHyGydZ/zaUZXtHkjOTvLyqrpfk/CTHJ7lrkvuveN8AAADAgC16RcfHkjy4qq6+fkE371dy5WteV6LNXhPzgCSvT/KHSc5I8hNJHtRaO32V+wYAAACGbdFGx0uS3C7J26vqjlV19W66Y2YPIb1tkhcvO+R6rbVvtdYe3Vq7cWvtmq21I1trb131fgEAAIBhW+jWldbam6rq2UmelORDmb3KtWXWMKkkz22tvWHpKQEAAADmsOgzOtJa+4OqemuShyW5VWYNjnOSvLa19pEl5wMAAACY28KNjiTpGhqaGgAAAMCgbKnRAQCrdMiug/uOAADASG3Y6Kiqp2b2DI5ntdau6D5vprXWnrGUdAAAAAAL2OyKjlMya3Q8N8ll3efNtCQaHQAAAMCO26zRcYskaa1dtvYzAAAAwBBt1ug4Psmb13xuSS5srX1vdZFg+06875F9RwAAAKAHB2yy/JQkP7Hm8/lJHriyNAAAAADbsFmjY3eSXWs+1wqzAAAAAGzLZreufCzJ46vqakn+tZt3t6racLvW2quWEQ4AAABgEZs1Oh6T2TM6Xth9bkl+q5v2pSXR6AAAAAB23GZXZny8qn4syS2T3CTJe5I8K8m7Vh8NAAAAYDEbNjqq6u5JPtta+1ySz1XVe5O8p7X23h1JBwAAALCAzR5GemaSe6z5fHiSQ1aWBgAAAGAbNmt0fD/JNdZ8vnmSQ1cXBwAAAGDrNnsY6blJjq+qs3LlW1euV1X/fqONWmtfXkY4gLG56Iw/6zsCAADs1zZrdDwzyWuTnNV9bkn+uJs2cuA2c7GgW936xn1HAAAAgN5t9taVv6qqjyc5KrO3rpyc5K1JPrH6aAAAAACL2eyKjux540qSVNUpSd7UWnvtinMBAAAALGzTRsdarbXNHl4KAAAA0JuFGh17VNXdk9wzyY2SvKC1dnZVHZrkyCSfaK3tXmJGYD/1tVN/u+8IALBfOmTXwX1HANiyha7QqKoDq+oNSc5M8uQkv57kpt3iH2T2/I5HLTUhAAAAwJwWvRXlCUmOTfKYJLdJUnsWtNYuTfKWJL+wtHQAAAAAC1i00XFckle11v4kyTf2svyzSf7DtlMBAAAAbMGijY7Dk/zTBst3J7nOltMAAAAAbMOiDyP9dpLrbrD8Vkku3HocAGAnHXP3I/qOAACwVIte0fGBJA+rqlq/oKquk9nDSc9cRjAAAACARS3a6HhWklsneXeS+3XzfrKqfivJWUkOSfKc5cUDAAAAmN9Ct6601j5aVQ9K8vIkr+hmPz+zt69ckOSBrbXPLDciAAAAwHwWfUZHWmvvqKrDk9wzyRGZNTk+l+SdrbVLlpoOAAAAYAELNzqSpLX2/SSndxMAAADAIGyp0VFV105ydJJbdrPOS/L3rbVvLysYAAAAwKIWbnRU1W8keUGSQzO7bSVJWpLvVNVjWmsvX2I+AAAAgLkt1Oioql9M8rLMruB4apJPdYtul+S/JXlZVV3QWnNLCwAAALDjFr2i4/FJPpvkTq2176yZ/w9V9YokH0zyhHh2BwAAANCDAxZc/yeTnLauyZEk6Z7P8cpuHQAAAIAdt2ijI7nyuRx707YaBAAAAGC7Fm10fDzJ8VV1yPoFVXVokhO6dQAAAAB23KLP6Hh+kjcnOauqXpTkM938PQ8jvVWSBy0vHgAAAMD8Fmp0tNbeWlWPTvLcJH+aK29VqSTfTfLo1trblhsRAAAAYD6LXtGR1tr/U1WvTXKPJLfIrMnxhSR/31q7eMn5AAAAAOa2cKMjSVpru5P85ZKzAAAAAGzLpg8jraoDq+o5VXXSJus9sqr+sKo2eisLAAAAwMrM89aVhyV5XJKPbLLeh5M8IcmvbjcUAAAAwFbM0+h4cJJ3tdb+50YrdcvfGY0OAAAAoCfzPKPjZ5K8YM7xzkzymK3HAQDoz2kn3avvCADANs3T6LhukgvmHO/Cbn0AAEbua6f+dt8RAGBh89y68u0k159zvOsl+c7W42yuqo6tqjdU1XlV9b2qOr+qXllVh69yvwAAAMDwzdPo+HSSe8453j269Vfp8UmumeTpSe6d5JQk/zHJWVV1ixXvGwAAABiweRodb05ydFXdf6OVquoXM2t0vGkZwTZwTGvt/q2101pr722tvTKzRsyuJI9e8b4BAACAAZun0fE/knw+yRur6lnrbxGpqsOr6plJ3pjk3G79lWmt/dDzQlpr5yf5RpKbrXLfAAAAwLBt+jDS1tr3quq+Sf4myZOSPLGqvp3kW0l+JMm1k1SSc5Lcr7V26Qrz7lVV3T7JDZJ8ah/Ld28yxGFLDwXAVajFAP1Sh4H9xTxXdKS19vkkP5Xkd5J8IMkPktw4yeVJ3t/NP7K19oUV5dynqrpGkpcn+WaSl+70/gEAAIDhmOf1skmS7kqNP+2mpaiqo5KcOefqN2itfWPd9gcmeVVmTZj7tdYu3NuGrbVdm+TYHR1sgJVSiwH6pQ4D+4u5Gx0rcnaSh8+57rfXfqiqA5K8IsmDkvxKa+3vl5wNAABgpY65+xF9R2Dgvnbqb/cdYXR6bXS01v45yWmLbtc1Of48yX9N8rDW2puXHA0AAADmcvbzju87Amv0fUXHwqqqkpya5NeSPLy19vqeIwEAAAADMbpGR5IXJfn1zJod51bVndcs+1Zr7TP9xAIAAAD6NsZGxzHd1xO7aa33JjlqR9MAAAAAgzG6Rkdr7fC+MwAAAADDdEDfAQAAAACWRaMDAAAAmAyNDgAAAGAyNDoAAACAydDoAAAAACZDowMAAACYDI0OAAAAYDI0OgAAAIDJ0OgAAAAAJkOjAwAAAJgMjQ4AAABgMjQ6AAAAgMnQ6AAAAAAmQ6MDAAAAmAyNDgAAAGAyNDoAAACAydDoAAAAACZDowMAAACYDI0OAAAAYDI0OgAAAIDJ0OgAAAAAJkOjAwAAAJgMjQ4AAABgMjQ6AAAAgMnQ6AAAAAAmQ6MDAAAAmIyD+g4AAAAASXLaSffqOwIT4IoOAAAAYDI0OgAAAIDJ0OgAAAAAJkOjAwAAAJgMjQ4AAABgMjQ6AAAAgMnQ6AAAAAAmQ6MDAAAAmAyNDgAAAGAyNDoAAACAydDoAAAAACZDowMAAACYDI0OAAAAYDI0OgAAAIDJ0OgAAAAAJkOjAwAAAJgMjQ4AAABgMjQ6AAAAgMk4qO8AAACwTBed8Wd9RwCgRxodAAAAwNL03XB26woAAAAwGaNvdFTVaVXVquqtfWcBAAAA+jXqRkdV3SPJLyX5Vt9ZAAAAgP6NttFRVYcmOTXJyUn+tec4AAAAwACMttGR5NlJLkryx30HAQAAAIZhlG9dqar/mOS3ktyltXZ5VfUdCQAAABiA0TU6quoaSV6e5MWttf855za7N1nlsG0HA2BDajFAv9RhYH/R660rVXVU98aUeabrd5udnORaSZ7SY3QAAABggKq11t/Oq26c5N5zrv66JLdM8okkxyd5x5pln0jyySQPTXJJa+2yBXNckaQOO0wTGxiWiy+++MuttZv3nWMnqMXAUO0vtVgdBoZq0Trca6NjUVX1gCRv2WS1R7bWXrrguD/I7OqWeV5Tu6fyX7zIPgZgjLnHmDkZZ+4xZk7GmXvRzBfvD39cJ/tFLR5j5mScuceYORln7jFmTtTivdoP6nAyztxjzJyMM/cYMyfjzL3SOjy2Rsf1k9x+L4ten+TzSf57knNba19bYYbdSdJa27WqfazCGHOPMXMyztxjzJyMM/cYMw/RGI/jGDMn48w9xszJOHOPMXMy3txDMtZjOMbcY8ycjDP3GDMn48y96syjehhpa+0bSd6zfn5VXZrkG621H1oGAAAA7D96fRgpAAAAwDKN6oqOfWmtHd53BgAAAKB/rugAAAAAJkOjAwAAAJgMjQ4AAABgMjQ6AAAAgMmo1lrfGQAAAACWwhUdAAAAwGRodAAAAACTodEBAAAATIZGxzpV9YqquqSqfmwvy15YVT+oqp/tPj+rqs6oqguqqlXVKTse+Mpsc+WuqutX1V9V1eer6jtV9a2q+lhVPbqqDhxi5u5z28d00hAzV9UJG2RuVfWQIebuPt+sql5TVd+squ9V1Uer6oEDyjfX711VHVNVr66qz1bV5VX1xZHkfnJ3zL9ZVd+vqvOr6mVV9aPLzj9kY6zF6vDwcg+pFg+9Dm8h4yBqsTq8Ourw8HJ3nwdRi8dYhxfJ3X32N3G/ubdXi1trpjVTksOSfCXJPyU5cM38uya5PMmz1sz7TpL/L8nLkrQkpww9d5KbJXl1kkckOTrJvZP8SZf/pUPM3M1rSV6f5M7rphsOMXOSG+wl652TfCjJ95LsGmju6yT5UpL/neSE7vz4iyRXJPmlvvN18+b6vUvy8iRnJ3lNknOSfHEkuf8oyROTHJPkqCSPTPLVJF9Lcp2dPG/6nFZxbIeSOerwTh7rwdTiBTL3Uoe3cD4MohavKLM6vKJjO5TMGVAd3sKxHkQtXuBYD6YOL5jb38T9595WLd6xk2pMU5L7dAf98d3nQ5J8Psknk1x9zXoHdF93bfRDGlrufWz7+iTfT3LQEDN36/xx3+fGdo5zkhsmuSzJa4aaO8mTuwL+0+u2PTPJl/ec8z2fC3P93q3NmuStqyjqq8i9j33cu9vmuD7Onb6mMdZidXh4ufeyXW+1eOh1eMHzYTC1WB0ez/kwpMz72LaXOrzgsR5MLR5jHZ43d5+1eIx1eBW597GPuWvxjp9YY5mS/HmSS5PcNsmLk/yfJEfuY93ei/pWcq/b7sWZdVVX+gfUVjMPqahv9Tgn+f3u+/j5oeZOcnqSr+xlu9/rst9lKMd1kd+7VRb1VeZes80du21+ta9zp69pjLVYHR5W7r1s02stHnodXvS4DqUWq8PjOx+Gknnddr3V4XlzD60Wj7EOz5O771o8xjq8ytxrtpm7FvdyYo1hyuzym/+d5AuZdfOetsG6gyjqi+ROUkkOyuyyrAdndgnR04eauTu+F3X/+Fya2eVuDx76cV63zaeTnJ+khpo7yTuTfGEv2z2q+xmcOJTjOrCivvTc3e/ntZL8VJIPJPlskkP6Onf6msZYi9XhYeXeyza91uKh1+FFj+tQarE6PL7zYQiZh1SH5809tFo8xjo8T+6+a/EY6/Cqcm+1FvdyYo1lSvIb3cE/J8nVlnFyDSV3kkd367TuJHzmkDNndh/lf01yt+4fovd06//OUDOvW/fO3bpPHfL5keSFSX6Q5N+tm/+abpsnDeW4DqmoLzt3kkPX/H62JB9McpO+z52+pjHWYnV4OLnXrTuIWjz0OrzIcR1SLVaHx3U+DCHz0OrwPLmHWIvHWIc3yz2EWjzGOrzs3Nupxd66sg9VdVBmHbsrktw8yQ89QXaIFsj9hiQ/m+SeSZ6T5Per6k93JOQ682RurT2stfba1tr7W2tvTPJzSd6f5JlVda0dDZwtnR+/3q172mqTbWyO3C/LrKi/tqpuV1XXq6r/luSXuuVX9JxvkFaQ+5LMfj//U2YPSrtOkvdU1U22Oe7ojPGcUId3zhhr8dDr8JwZB0cdXp2Jnw+DqcPJOGvxGOtwMvxaPMbfu2RgtXinO2djmZI8tfsB3T+zB858OGueILtoN2qIuddt99jue/jpVWdcYubf7DL/7JAzJzk4ycVJ/m4M50dmDxL6cq7snH45V/4fj1/rO9+adQfTvV51vUhy08wejvYnfZ9DOz2NsRarw8PMPZRaPPQ6vIXjOoharA6P73wYSuZ12/VWh7eZ29/EK8jdZy0eYx1eZe4128xdi3s9wYY6JfmJzJ4E/OLu872y5gmyy/ghDSH3um3v1q37KyPKfFK37s8MOXOS4/o4ttvJndk9q7dOcpskB2Z2ieQVSW4xhHzd8kEU9Z2qF5nd63hGn+fQTk9jrMXq8HBzD6EWD70Ob/G49l6L1eFxng9Dybxu217q8BJy+5t4Rbn7qMVjrMOrzr1uu7lqcW8n2FCnJFdL8rEk52XNQ04ye//w95IcsawfUt+5123/1J0ukNvJnOSAzC7T+1aSaw45c2b3Tn4zyTXGeH4kuUZm3di/GVK+IRT1naoXSW6Z2eWTf9rXObTT0xhrsTo87GPddy0eeh3exnHttRarw+M9H4aSed32O16Ht5u7r1o8xjq8hGPtb+Kecq/ZZu5a3MsJNuQpycmZden+y7r5e54g+4+58t2//zmz+7T2dCff2H3+pSQHDzF3Zq9zenmShyY5KrPLil7SnTB/NeDMpyb51S7zQ5K8tzvmjxpi5jXzb9mt3+sfRgsc6wOSvCjJg7pj/Ygkn8js0rMf7TtfN2+u37vM7gvcM//DSS5Y8/m2Q8zdZX5/ZpdF3jvJ0Ul+tzv+FyY5vM/zaIjn7CLnxFAyRx3e0fOjm997LV7gWPdShxc9rvP+3mXFtXjZmaMOr/R8GErmDKgObyH3IGrxIudHN7/3OrzgsfY3cY+5s4Ra3NtJNsQpyU9mdrnNS/ax/D7dD+Ox3ef35KpPgV07bXrw+8jdnSR/m+Tr3TbfTvKRJL+T5KCBZj6mO9EvzOxdzLuT/EOSY4Z8fnTznpEe7/XcwrE+IMlfJ/lat81Xkrw0yY2HkK/7PNfvXZITNljvlCHmzuwfg9OSnJvZK+6+n1ln/H8kuXlf59CQz9lFzomhZI46vKPnRzev11q84LHe8Tq8leM67+9dVliLV5E56vBKz4ehZM5A6vAWcg+iFi96fnTz/E285OM67+9d9tO/iasbCAAAAGD0vF4WAAAAmAyNDgAAAGAyNDoAAACAydDoAAAAACZDowMAAACYDI0OAAAAYDI0OgAAAIDJ0OhgFKrqqKpq66bvVNVZVfV7VXVQ3xkBpk4tBuiXOgzz8YvA2LwuyTuSVJIbJzkuyf+d5DZJfrPHXAD7E7UYoF/qMGygWmt9Z4BNVdVRSc5M8rjW2vPXzD8kydlJ/l2SG7XWLuwh24FJrtFau2Sn970sVfUjrbVv950DGDa1eLXUYmAz6vBqqcPT4dYVRq219t0kH8ysm/0f1i6rqjtW1Vuq6htV9f2qOqeq/mBvl/RV1bFV9fGqurSqvlxVJ1fV0d3lgCesWe+Ebt7RVfWUqvpCkkuTPHjR/VbV7arqL6vqq916/1xVZ1bVfdesc82qOqUb45Kq2l1Vn6yq5+3le/iN7rLF71XVxVX1d1V1172s16rqtKr6+ar6QFV9J8npCxx2gKtQi68ynloM7Dh1+CrjqcO4dYVJ2FPML9ozo6p+Iclbknw+yQu6ZXdJ8vQkP5Xkl9es+yuZXf73hSRPS/KDJMcnOWaDfT4/ydWSnJrkW0nOWWS/VXW9JO/uxnppki8luX6SOya5U5K3d8tekuTXk7wqyQuTHJjk1kl+bm2Yqnpukscn+XCSJyf5kcwuWzyzqu7fWnvHuvx3THJsl/+VG3yfAPNSi9VioF/qsDrMHq01k2nwU5KjkrQkT82s+N0gyR0yK3otyYfXrHvNJP+c5H1JDlo3zu916x/VfT4oyVeT/EuS66xZ79Ak53XrnrBm/gndvHOSHLxu7EX2+4vd5wdv8n1flOQdm6zz40muSPKBJFdfM/+mSXYn+WKSA9fMb910dN8/V5PJNK5JLd5wHbXYZDKtfFKHN1xHHTb92+TWFcbmaUkuTHJBkk8keVSSN2dWJPe4R5IbJXlFkl1Vdf09U2YPbUqSe3Zffyaz4ndaa+1f9wzQWvtOZl3lffmz9sP3Hy6y34u7r/epqmtvsJ+Lk9yuqm6/wTr3z+wyxT9qrV225nv4WpLTktw8yU+v2+bjrbV3bTAmwEbU4h+mFgM7SR3+Yeow/8atK4zNy5L8ZWaXyN0hyROS3CyzewL3uE339c83GOdG3ddbdF/P2cs6e5u3x7l7mTf3fltr762qV2XWDX9oVX0kybuSvKG19pk16/9ukr9I8smqOi+zh0+dnuT01toV676HT+9lf5/qvt4yyUc3yQ8wL7VYLQb6pQ6rw2xAo4Ox+dyarusZVfWBzC5Pe2mSh3Tzq/v6uCT/ax/jfG3duova29OkF9lvWmvHdw9Q+oUkd03y2CR/UFW/21p7cbfO26rq8G6d/5zk6CSPSPL+qjq661Zv5XsY7dOwgUFQi9VioF/qsDrMBjQ6GLXW2j9W1V8kOa6qXtRa+8ckn+sWf3eOS9HO777++F6W7W3eRhbZb5KktfapzDrMf1RVu5J8KMlzquolrXU3D7Z2UZJXJ3l1VVWS52T2kKX7Z9bJ/0I33O3W/Pcet+2+nrfg9wIwN7VYLQb6pQ6rw1yVZ3QwBc9IcnlmT3FOkndmdr/iE6vquutXrqprVdWPdB8/muTrSU6oquusWefQJCctmGPu/VbVdavqKr9/rbXdmf0jc3CSa1bVgV2hX7tOS/Kx7uOeffx1Zg9SelxVXW3N/m6S5OGZPb36YwFYLbVYLQb6pQ6rw3Rc0cHotdY+X1Wvz+y+vru11t5fVccleWuSc6rqzzN7tdWuJEckeVCSByZ5T2vtB1X1+0lek+TDVfXyzF6ldUKSb2Z2r1+bM8d3591vkuOS/F5V7Xnt1v/J7DK8eyV5Y2vte11B/3pV/XVmRfmCLs8jk/xruvd8t9bO6S73e3yS91XVG3Llq7QOTfLQ1trlix5XgEWoxWqZUPvJAAABY0lEQVQx0C91WB1mjb5f+2IyzTPlyldp/f4+lt8msw72mWvm3T6zy9u+muSyzF6X9Y9JnpLkuuu2f3BmT6z+fpIvJzk5swJ8lddd5cpXaR21QdZN95vZ+8NfmVlB/25m7x3/eGb3JF6jW+fqSZ6d2XvAv9ll+2JmD3a69V72e2Jmxf/Sbry/T3K3vazXMnuidu8/V5PJNK5JLVaLTSZTv5M6rA6b5puq+yED61TVY5M8P8ldWmsf7DsPwP5ILQbolzrMGGl0sN+rqqsnubytuZStux/xE0muneSmbc27uAFYPrUYoF/qMFPiGR0we5/2Gd09jecnuUmS49Pd+6egA+wItRigX+owk6HRAcmFST6Y5KFJbpjZg5c+meSJrbU39hkMYD+iFgP0Sx1mMty6AgAAAEzGAZuvAgAAADAOGh0AAADAZGh0AAAAAJOh0QEAAABMhkYHAAAAMBkaHQAAAMBk/P/3FHQO674WUgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1296x720 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Xs = ['X'+str(i + 1) for i in range(X.shape[1])]\n",
    "lams = [10**4, 10**2, 0]\n",
    "\n",
    "fig, ax = plt.subplots(nrows = 2, ncols = len(lams), figsize = (6*len(lams), 10), sharey = True)\n",
    "for i, lam in enumerate(lams):\n",
    "    \n",
    "    ridge_model = RegularizedRegression()\n",
    "    ridge_model.fit_lasso(X, y, lam) \n",
    "    ridge_betas = ridge_model.beta_hats[1:]\n",
    "    sns.barplot(Xs, ridge_betas, ax = ax[0, i], palette = 'PuBu')\n",
    "    ax[0, i].set(xlabel = 'Regressor', title = fr'Ridge Coefficients with $\\lambda = $ {lam}')\n",
    "    ax[0, i].set(xticks = np.arange(0, len(Xs), 2), xticklabels = Xs[::2])\n",
    "    \n",
    "    lasso_model = RegularizedRegression()\n",
    "    lasso_model.fit_lasso(X, y, lam) \n",
    "    lasso_betas = lasso_model.beta_hats[1:]\n",
    "    sns.barplot(Xs, lasso_betas, ax = ax[1, i], palette = 'PuBu')\n",
    "    ax[1, i].set(xlabel = 'Regressor', title = fr'Lasso Coefficients with $\\lambda = $ {lam}')\n",
    "    ax[1, i].set(xticks = np.arange(0, len(Xs), 2), xticklabels = Xs[::2])\n",
    "\n",
    "ax[0,0].set(ylabel = 'Coefficient')\n",
    "ax[1,0].set(ylabel = 'Coefficient')\n",
    "plt.subplots_adjust(wspace = 0.2, hspace = 0.4)\n",
    "sns.despine()\n",
    "sns.set_context('talk');"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
